Commodity Asian option pricing and simulation in a 4-factor model with jump clusters

Mean reversion, stochastic volatility, convenience yield and presence of jump clustering are well documented salient features of commodity markets, where Asian options are very popular. We propose a model which takes into account all these stylized features. We first state our model under the historical measure, then, after introducing a structure preserving change of measure, we provide a risk-neutral version of the same model and we show how to price geometric and arithmetic Asian options. To this end, we derive semi-closed formulas for the geometric Asian options price and develop a computationally efficient simulation scheme for the price process, allowing to price the arithmetic counterparts using control variate technique. Finally, we propose a simple econometric experiment to document presence of jump clusters in commodity prices and evaluate the performances of the proposed simulation scheme on some parameter sets calibrated on real data.


Introduction
Commodity derivatives markets had a tremendous growth in recent years, both in trading volume and variety of offered products, raising the need for models able to replicate correctly observed market prices. To quantify this growth, we present in Fig. 1 the time series of trading volumes of commodity futures and options. The total number of trades has more than doubled over the last ten years. Moreover, we note that in 2020 volumes increased substantially with respect to 2019 (by 35.7% for futures and 26.2% for options). This is possibly due to the fact that during highly uncertain times (like the COVID-19 outbreak) investors are more likely to hedge risk than in calm periods (see Gonzato and Sgarra, 2021 and references therein).
Much literature has been devoted to the study of empirical properties of commodity prices. Bessembinder et al. (1995) find clear evidence of mean reversion across many commodity markets. This is also confirmed by other prominent studies such as Schwartz (1997), Casassus and Collin-Dufresne (2005), and others. In commodity markets mean reversion is mainly induced by convenience yields, which stem from both the reduction in cost of acquiring inventory and the value of being able to profit from temporary local shortage of the commodity (Yan, 2002). Lutz (2010) surveys different methodologies proposed in literature to jointly model mean reversion and convenience yield and concludes that the so-called "autonomous convenience yield" approach provides good empirical performances for a wide variety of commodities. Based on this finding, we model mean reversion in the spot price directly by means of an Ornstein-Uhlenbeck process and assume that the convenience yield follows a mean reverting process which is independent of the spot price. Since typically convenience yield may assume both positive and negative values, following Casassus and Collin-Dufresne (2005), we assume a Gaussian Ornstein-Uhlenbeck process. Another salient feature of commodity markets is stochastic volatility of the log-returns, as documented, among others, by Trolle and Schwartz (2009) and Cortazar et al. (2017). Finally, there is intuition for the presence of self-excitation in the jumps process. This is the phenomenon, also named "jump clustering", in which whenever a jump in the asset price occurs, the possibility of observing subsequent jumps increases. Filimonov et al. (2014) found clear evidence of this effect in  (2019) for data until 2019 and https://focus.worldexchanges.org/articles/commodity-derivatives for data on 2020 many commodity markets and concluded that: "At least 60-70 per cent of commodity price changes are now due to self-generated activities rather than novel information". Additional studies confirm the importance of including self-excitation in commodity price modeling for derivatives pricing, hedging and forecasting (Jiao et al., 2019 andGonzato and. In the course of this paper we provide further evidence of jump clustering effects in oil and precious metal markets. Following standard asset pricing literature (see Fulop and Li, 2019 and the references therein), we model the stochastic jump intensity through a Hawkes process with exponential kernel. The resulting proposed model is a 4-factor affine model for commodity prices. The model takes into account all aforementioned stylized features and its primary scope is to accurately price commodity derivatives. To this end, after describing our model under the historical measure, we propose a structure preserving change of measure which allows formulation of the model under a risk-neutral pricing measure. The necessity of using general multifactor models for pricing commodity derivatives is highlighted in Cortazar et al. (2017) and Schöne and Spinler (2017). The model can be considered as an extension of that proposed in Casassus and Collin-Dufresne (2005) with stochastic volatility and self-exciting jumps. Due to its affine structure we obtain closed formulas for the futures price and semi-closed formulas for European options on futures, useful for model calibration purposes, as detailed in Sect. 5. Moreover, as we discuss later on, the process can be accurately simulated and allows for efficient evaluation procedures for Asian options. This is very important for several reasons: (i) in commodity markets Asian options are traded much more frequently than in equity or interest rate markets; (ii) model sophistication typically precludes exact formulas for the price of exotic derivatives and pricing must be performed through simulation; (iii) general multi-factor models present multiple sources of randomness leading to a large variance in the simulation step and, if the simulation is not accurate enough, to relevant mispricing of the derivative instrument (as we document in the case of the Euler scheme).
Asian options are very popular among commodity derivatives traders and risk managers,Here we don't mean that Asian options are the most traded options in general, but that they are strongly linked to the commodity market and they are quite popular among practitioners, especially in OTC markets (see e.g. Kaminski, 1999). as their payoff depends on the arithmetic average of the prices assumed by the underlying during the life of the contract and represent a cheaper alternative to European options. The averaging process is responsible for their popularity: it smooths possible market manipulations occurring near the expiry date, reduces payoff's volatility and allows for better cash flow matching. Hence, many institutions started quoting average price options for highly volatile assets such as commodities (we refer to Roncoroni et al., 2015, Chapter 18 for a nice description of the traded products). However, Asian options are traded mainly in Over the Counter (OTC) markets (e.g. Kaminski, 1999).
Even if exact computation of the price of an arithmetic Asian option is precluded, we obtain semi-closed formulas for the price of a discretely monitored geometric Asian option under the proposed model. In this context, we contribute to the literature deriving the moment generating function of the arithmetic average of the asset returns, extending the results in Fusai and Kyriakou (2016) to the case of multi-factor affine models with mean reversion. Then, for accurate model simulation, we rely on two different streams of literature on the exact simulation of option pricing models with stochastic volatility (Broadie and Kaya, 2006) and on the exact simulation of affine point processes (Dassios and Zhao, 2013). We use their results as building block for the construction of an almost exact simulation scheme for the proposed model which allows to simulate asset price trajectories observed at discrete times and produce almost unbiased estimates for exotic derivatives prices. Unfortunately, exact simulation is possible only in absence of mean reversion, in the opposite case we propose a simple approximation, whose accuracy is investigated through extensive numerical experiments. Due to the approximation, accuracy increases with the number of time discretization steps. Since we consider four stochastic factors we observe high Monte Carlo variance in the simulation step. Therefore, we employ the previously computed geometric Asian option price as control variate to reduce the variance of the estimator for the arithmetic Asian option price (Kemna and Vorst, 1990). Effectiveness of the proposed pricing methodologies is confirmed through extensive numerical experiments using realistic model parameters calibrated on real market quotes. We show that the bias of the proposed simulation scheme decays faster than that of the Euler scheme (used as benchmark) with a higher convergence rate of the Root Mean Squared Error (RMSE) and better overall performances in terms of trade-off between accuracy and computational effort.
The rest of the paper is organized as follows. In Sect. 2 we introduce the model specification and derive pricing formulas for European call options on futures (which will be later used for model calibration). In Sect. 3 we present an efficient simulation algorithm for the proposed model, while Sect. 4 deals with Asian option pricing. In Sect. 5 we perform numerical studies, where we document the presence of jump clusters and investigate the performances of the proposed simulation scheme. Section 6 concludes.

Model setup
In this section we outline our model under the historical probability measure. Then, we introduce a structure preserving change of measure and describe our model under the riskneutral measure. Finally, we show how to price futures and options on futures contracts under this model specification.

Model dynamics under the historical measure
Let ( , F , (F t ) t∈ [0,T ] , P) be a filtered probability space, which supports all the processes we encounter in the sequel. If we denote by S t the spot price process, the log-return process is X t := log(S t /S 0 ). Under the historical measure P the dynamics of the log-returns is defined by the following system of stochastic differential equations: where J x N t is a marked point process with stochastic intensity λ t and random jump size J x and μ * = e μ J +σ 2 J /2 − 1 its compensator. We assume a set of initial conditions satisfied by each one of the processes described by the SDE system: {X 0 , V 0 , δ 0 , λ 0 }. Following Casassus and Collin-Dufresne (2005, Proposition 2), we assume that the convenience yield δ t evolves according to an independent Gaussian Ornstein-Uhlenbeck process (which can assume both positive and negative values). The parameter α ≥ 0 controls the speed of mean reversion. Following Larsson andNossman (2011), Brooks andProkopczuk (2013), Gonzato and Sgarra (2021), we assume a Cox et al. (1985) (CIR) square-root diffusion process for the variance V t , which is correlated with the price process through the coefficient ρ ∈ [−1, 1], i.e. E dW x t dW v t = ρdt. We model abrupt changes in log-returns X t by including an independent compound Poisson process with stochastic intensity λ t , such that the counting process exhibits a self-exciting behavior, which means that a jump increases the probability of observing subsequent jumps. Although the jump intensity λ t is stochastic and path-dependent, it is possible to prove that the vector process {X t , V t , δ t , λ t } is Markovian and affine (as we discuss later on). The characterization of price jumps is completed by specifying a probability density function for their sizes, which are i.i.d. Gaussian and denoted by J x ∼ N (μ J , σ 2 J ). The coefficient μ describes a constant drift with respect to the historical measure. We further assume that the following non-explosion condition for the Hawkes process N t in (4): and the well-known Feller condition on the parameters of (2) (granting the strict positivity of V t for all t > 0) 2k v θ v ≥ σ 2 v hold. The model nests many other models proposed in literature to describe a large variety of commodities dynamics, we mention just a few of them. The classical Gibson and Schwartz (1990) model is obtained with V t = λ t = 0 and constant interest rates, the model in Casassus and Collin-Dufresne (2005), well suited for a large number of commodities such as oil, copper, gold, silver is obtained imposing V t = λ t = 0. Eydeland and Geman (1998) propose a model for gas and electricity which is obtained with δ t = λ t = 0, while the one of Geman (2000) (for oil) with λ t = 0. Moreover, we mention the model proposed in Larsson and Nossman (2011) for oil prices which is a special case with δ t = α = 0 and constant λ t .
Before performing our analysis on the proposed model we want to examine the issue related to existence and uniqueness for the system of stochastic differential equations characterizing our model. We have the following Proposition 1 There exists a solution to the system of (1-4) when a non explosion condition holds and this solution is adapted to the filtration generated by the driving processes. This solution is unique when standard initial data (F 0 -measurable initial data) are assigned.
Proof See Appendix A.
As far as positivity of (2) and (4) is concerned (neither positivity of X t nor positivity of δ t are required at all), the classical Feller condition on the parameters describing the volatility dynamics 2k v θ v > σ 2 v will grant the strict positivity of V t for t > 0, while (4), as an Ornstein-Uhlenbeck SDE driven by a jump process with positive long-term mean θ λ and with only positive jumps in the driver, cannot exhibit a negative solution.

Risk-neutral dynamics
In order to deal with derivatives pricing we need to introduce a risk-neutral measure. To this end, we extend the results in Zhang et al. (2009) and Hainaut and Moraux (2018) by introducing an Esscher-type measure change. This is defined as follows. If we denote by L t := J x d N t the jump term in (1), and by ψ P (z) := E[e z J x ] the moment generating function of the jump size density, we can define the following family of exponential martingales: where κ 1 (ξ ), κ 2 (ξ ), η (with κ 1 (ξ ), κ 2 (ξ ) functions of ξ ) denote the risk premium related to the jumps and φ x , φ v , φ δ (stochastic processes adapted to the reference filtration F t ) denote the risk premium related to the diffusion components of log-returns, volatility and convenience yield respectively. We state the following result.
Proposition 2 If, for any ξ , there exist functions κ 1 (ξ ), κ 2 (ξ ), which are solutions of the following system of algebraic equations: is a local martingale. If, moreover, the non-explosion condition (5) holds for N t and the Novikov condition holds for φ x , φ v , φ δ , M t is a true martingale.
Proof See Appendix B The following result shows that the measure change introduced by the likelihood process preserves the model structure.

Proposition 3 The dynamics under the risk-neutral measure
δ Q t is given by the following system of stochastic differential equations: where J Q x N Q t denotes the jump process with respect to Q and μ ,Q = e μ Q J + 1 2 σ Q,2 J − 1. The relations between the relevant parameters under P and Q are the following: Moreover, under Q the jump size J Q x is still normally distributed with mean μ Q J and volatility σ Q J , with moment generating function given by ψ Q (z) = ψ(z + ξ)/ψ(ξ). Proof See Appendix C.

Remark 1
We shall assume in the following that the risk premium terms φ x , φ v are such that the following equality is satisfied for all t ∈ [0, T ]: in such a way that the dynamics of log-returns under Q can be written as in (7), by observing that ψ(1) − 1 = μ * .

Remark 2
The last sentence in Proposition 3 implies that the relations between the mean and the variance of the jump size with respect to Q and P are the following: We point out that when α > 0 the spot price does not satisfy the standard no-arbitrage condition and the spot price process is not a local martingale with respect to Q. This is not a problem since commodities with storage costs of the good are not directly traded assets, therefore the drift of the log-returns can be of the mean-reverting type under the risk-neutral measure (see e.g. Schwartz, 1997;Lutz, 2010;Benth, 2011;Cai et al., 2014). As a result, the proposed framework is well suited also for option pricing, as witnessed also by the large amount of literature dealing with option pricing under mean reversion of the asset price (see Fusai et al., 2008;Wong and Lo, 2009;Chung and Wong, 2014;Brignone et al., 2021, among others). However, in the special case with α = 0 the price process is a local martingale under the risk-neutral measure Q. We are also assuming that the short rate is equal to 0. This is in agreement with most theoretical models with stochastic convenience yields (e.g. Schwartz, 1997;Routledge et al., 2000).

Remark 3 The non-explosion condition under Q is now given by
that we'll assume to be satisfied in the following.
In the next subsection, we show how to perform option pricing under the model (7)-(10). Despite the growing interest and success in modeling jump clustering in finance (see e.g. Fulop and Li, 2019 and the references therein), the application of such framework on commodity options is still unexplored and we provide a first contribution on this topic.

Remark 4
Since in the rest of the paper we shall work always under the risk-neutral dynamics, we shall drop the superscript Q from all the relevant quantities.

Pricing options on futures contracts
We derive, next, the joint moment generating function of the quantities described in (7)-(10), we will make use of this result to price derivatives under the proposed model:

Proposition 4 Given a final date T > t and the time to maturity
Moreover, in the case with α = 0, we have Proof See Appendix D.
Remark 5 C and G can be computed full explicitly and the solution is reported in Proposition 4. In addition, also B can be computed explicitly, however, when α > 0 the solution is in terms of hypergeometric functions. Hence, we find in practice more convenient to solve the corresponding Ordinary Differential Equation (ODE) numerically than through the analytic solution. D is discussed in the next remark.

Remark 6
The couple (N t , λ t ) is an affine Markov process (see e.g. Errais et al., 2010). Therefore, the model (7)-(10) is affine and it is possible to derive an expression for the moment generating function of log-returns as solution to a ODEs system. Nevertheless, an explicit solution for the ODE F 4 in Proposition 4 is not available and numerical solvers must be considered (see e.g. Errais et al., 2010;Da Fonseca and Zaatour, 2014).
In the following, we will exploit Proposition 4 to perform derivatives pricing. Consider S T = S t e X T where τ = T − t is the time to maturity and S t is the price of an asset at time t. The price of futures contract F(t, T ) is given by the standard relation: Under the proposed model, we have where we compute the expectation by replacing u 1 = 1 and u 2 = u 3 = u 4 = 0 into (11).
We turn now our attention to the pricing of European options. Consider a maturity T and a strike K , then the price will be given by where f X T (x) is the probability density function of X T obtained by numerical inversion of the characteristic function of X T , i.e. E[e iu X T |F t ], a special case of Proposition 4 with u 1 := iu and u 2 = u 3 = u 4 = 0. Given a (univariate) characteristic function we compute the corresponding probability density function through the Fourier-Cosine (COS) method proposed by Fang and Oosterlee (2008). Given the probability density function, the European option price can be obtained in several ways, for example by solving the remaining integral (e.g. 13) using the trapezium rule.
Consider now a European option on futures at the initial date t, with maturity of the option T and maturity of the underlying futures contractT > T , the price is given by: where

Proposition 5 The moment generating function of Y T is given by
Given the moment generating function in (15) we get the characteristic function by substituting u with iu, then we compute f Y (y) using the COS method and recover C E F from (14). In the next section we tackle the problem of efficiently simulating the model in (7)-(10), while the pricing of Asian options is deferred to Section 4. X t j for j = 0, . . . , n. Anyway, for sake of clarity, let us illustrate the proposed method in the case where one wants to simulate directly (X T |X 0 , V 0 , δ 0 , λ 0 ), omitting intermediate dates. This is without loss of generality since only simple adaptations are needed to include also the intermediate dates, as we show in Algorithm 1 where we summarize the whole simulation procedure for a generic set of dates. Given the underlying asset price trajectory is then possible to price many kind of exotic options on spot prices, including Asian options. For what concerns the price of options on futures, we note that the underlying can be simulated easily from (12) requiring only the additional simulation of {δ t j } n j=1 (not required in the case of the spot price).
Step 1: Exact simulation of T 0 s ds given 0 Using Dassios and Zhao (2013, Algorithm 3.1) we obtain a sample of the triplet where N T is the total number of jumps in the period [0, T ] and τ k is the k−th jump time. Given the triplet we can compute Step 2: Exact simulation of ı T , T 0 ı s ds given ı 0 where The transition of δ T , T 0 δ s ds can thus be simulated exactly and efficiently by means of standard random numbers generators from a multivariate normal distribution.
Step 3: Exact simulation of V T , T 0 V s ds given V 0 The transition density of the terminal variance is known where χ 2 d (λ) denotes the non-central chi-squared distribution with d := 4θ v k v /σ 2 v degrees of freedom and non-centrality parameter λ : is simulated exactly using standard generators from the non-central chi-squared distribution. 1 The next step consists in simulating T 0 V s ds|V T , V 0 , not a trivial problem. Broadie and Kaya (2006) develop an exact simulation scheme. Nevertheless, their proposed methodology presents several implementation issues, being generally slow to run and, despite theoretically exact, biased in practice (see e.g. Glasserman and Kim, 2011;Kienitz and Wetterau, 2012). To overcome these issues, Kyriakou et al. (Forthcoming) suggest an alternative approach based on fast computation of the moments of T 0 V s ds|V T , V 0 and subsequent random sampling from a 4-moments matched Pearson distribution. This methodology turns out to be faster and more accurate than competing benchmarks, hence, we will adopt this methodology for fast sampling Step 4 From (7) and (19): is a sequence of i.i.d. normal (with mean μ J and standard deviation σ j ) random variables. Thus, we have wherē Therefore, we can simulate {J x,i } N T i=1 and the quantity X T + α T 0 X s ds by means of standard random numbers generators from a normal distribution. Note that in case of absence of mean reversion (α = 0) the model is simulated exactly through the proposed approach. However, in presence of mean reversion, a full exact simulation scheme seems not doable. Hence, we employ a central discretization (compare with Andersen, 2008, Eq. 32): If we denote with Z the random sample from X T + α T 0 X s ds, then we have that X T ≈ Z −αT /2X 0 1+αT /2 =:X T .
The whole simulation procedure contains two sources of error: i) the simulation of T 0 V s ds|V T , V 0 ; ii) the approximation in (21). The first is negligible in practice as illustrated in Kyriakou et al. (Forthcoming) and also confirmed by the numerical studies we will present in Sect. 5. The second source of error is more important. In particular, accuracy depends on the length of the interval [0, T ] (we expect high accuracy for small values of T ) and on the value of the parameter α. In this case, the smaller α the smaller the error. We will investigate this point in the numerical section.
Finally, the asset price is computed as S T = S 0 e X T . In addition, it is possible to recover T 0 X s ds from (21) and the price of futures from (12).

Asian option pricing
In this section we derive formulas for the price of Asian options. The payoff depends on the average of a commodity's spot-price. Consider the usual set of dates 0 =: t 0 < t 1 < · · · < t n := T , then the payoff of an Asian option with strike K and maturity T is where the arithmetic and geometric averages are defined respectively as A n = 1 n n j=1 S t j and G n = exp 1 n n j=1 log S t j . This kind of option is very popular in OTC markets, especially when the underlying is the price of metal commodities (see, among others, Kaminski, 1999, Shiraya andTakahashi, 2011). In practice averaging is usually arithmetic rather than geometric. However, pricing geometric Asian options is still very important as i) this allows for exact pricing formulas and ii) the payoff of geometric averaged options is highly correlated with that of their arithmetic averaged cousins. This fact is exploited in synthetic variance reduction methods for Monte Carlo simulation. Consequently, it comes natural to use the geometric Asian option price as a control variable in a Monte Carlo simulation to obtain accurate price estimates for the arithmetic counterpart. Let's start our discussion from the price of the geometric Asian option which is given by: where H T := 1 n n j=1 X t j . We derive, next, the moment generating function of H T . The following proposition extends the results in Fusai and Kyriakou (2016) to the case of a mean reverting multi-factor affine model.
Proof See Appendix E.
Given the moment generating function we obtain the characteristic function replacing u with iu and we price the geometric Asian option using the COS method. We turn now our attention to the pricing of the arithmetic counterpart. Kemna and Vorst (1990) first proposed to exploit the high correlation between the arithmetic and the geometric average of asset prices to reduce the variance of the Monte Carlo simulation estimator. In particular, the price of the geometric counterpart is used as control variable. We employ in our context a standard Control Variate Monte Carlo setup where, for a high accuracy, we simulate the underlying asset price process using Algorithm 1. We summarize the whole procedure in Algorithm 2 and refer to Glasserman (2004) for a detailed description of the usage of control variates method for Asian option pricing.

3:
Simulate {X t j } n j=1 using Algorithm 1 4: (22) 6: end Simulate {X t j } n j=1 using Algorithm 1 10: Compute S t j = S t 0 e X t j for j = 1, . . . , n 11: Compute X i and Y i from (22) 12: end We have considered options written on the average of spot prices. Another possibility is to average futures rather than spot prices. For sake of brevity, we will not consider this case in this work. However, the payoffs of the Asian options on futures can be obtained replacing S t j with F(t j ,T ) into (22). Results obtained for the spot price can be extended to the case of Asian options on futures similarly to the case of the European options (see Proposition 5). As a final remark, we point out that is also possible to price continuously monitored geometric Asian options under the proposed model applying the method outlined in Hubalek et al. (2017, Proposition 3) and extended in Brignone and Sgarra (2020, Proposition 2).

Numerical results
In this section we present numerical results. Before assessing the accuracy of Algorithm 1 and the pricing of Asian options we start by examining the time series of spot returns for four different commodities, namely, WTI crude oil, gold, silver and copper. We detect jumps in the historical price time series and present an econometric test which shows that jumps are not uniformly distributed over time, but tend to appear in clusters, providing further support for the proposed model. Then, we calibrate the model on real market option prices. Finally, using calibrated parameters, we evaluate the accuracy of Algorithm 1 and present numerical results on the pricing of Asian options.
All the computations are done using Matlab ® (Version R2021a) in Microsoft Windows 10 ® running on a machine equipped with Intel(R) Core(TM) i7-9750HQ CPU @2.60GHz and 16 GB of RAM.

Jump clusters in commodity prices
We illustrate and discuss the clustering effects of commodity prices, which provides us with an empirical justification for the introduction of self-excitation in the jump process. To this aim, we consider the historical time series of the spot prices from 30-Aug-2000 to 11-Dec-2020 (5048 daily observations) of four different commodities: gold, silver, crude oil, copper. 3 The sample involves periods of crisis and financial turmoil such as the 2008 credit crisis, the summer 2011 European sovereign debt crisis, and the recent COVID-19 outbreak. Spot prices are displayed in Fig. 2 along with jump occurrences (black bars). In order to detect jumps we employ the iterated re-weighted least squares technique developed in Callegaro et al. (2017) (see also Bernis et al, 2021 for more implementation details). Jump occurrences are displayed in Fig. 2. We note that jumps are very frequent in commodity markets: we identify in the whole sample a total of 144 jumps for gold price, 210 for silver, 137 for crude oil, 135 for copper. Therefore, a model which omits the jumps in the price process is likely misspecified. We also observe that jumps appear in clusters. This is particularly evident in the case of crude oil, where we observe prolonged periods of tranquillity (e.g. we observe no jumps between 29 Jun 2012 and 26 Nov 2014) followed by periods with a lot of jumps (11 jumps in 2015). In order to show that jumps do not exhibit a constant arrival rate we propose a statistical test. If the arrival rate is constant then the jump process is a homogeneous Poisson process and the distribution of the interarrival times is exponential with mean 1/λ, where λ is the average arrival rate computed as the ratio between the total number of jumps and the number of observations in the sample (e.g. 144/5048 = 0.0285 for gold). We then perform a two sample Kolmogoroff-Smirnov (KS) test, where null hypothesis is that the interarrival times are exponentially distributed with mean 1/λ. The null hypothesis is rejected at the 5% significance level for all the commodities considered, supporting the idea of stochastic jump intensity. In particular, for crude oil we obtain a p-value of 9.30E-07 which indicates a strong rejection of the hypothesis of a constant arrival rate. A graphical illustration is provided in Fig. 3 where we compare the empirical distribution of the interarrival times with the theoretical exponential distribution. If the arrival rate was constant then the two cdfs would match, but we note that this is not the case. The choice of using an Hawkes process to model the jump intensity (which is standard in the literature on equity price modeling) allows to take into account this feature of commodity returns.

Calibration
We calibrate the proposed model (7)-(10) on observed market option prices. In commodity markets the most liquid quoted options are American options. Hence, for a given date (Friday 5 March 2021) we collect American option prices for four different commodities: crude  Empirical (blue points) and theoretical (red lines) cumulative distribution function (cdf) of the interarrival times for four different commodities. The theoretical cdf is the one of an exponential distribution with mean 1/λ, where λ is the daily average arrival rate (0.0285 for gold, 0.0416 for silver, 0.0271 for crude oil, 0.0303 for copper) oil WTI, gold, silver and high grade copper. 4 We consider three different maturities for each commodity ranging from 1 to 9 months (options with longer maturities are not liquid enough). In order to avoid liquidity problems, we keep all the contracts whose trading volume at the end of the day is greater than two. We end up with 125 options for crude oil WTI, 57 for gold, 44 for silver and 38 for copper. Then, we convert American option prices to Europeans. Following Trolle and Schwartz (2009) we compute the implied volatility from American option prices using the Barone-Adesi and Whaley (1987) formula, then, we calculate European call options on futures prices using the Black (1976) formula. Hence, we can calibrate the model by matching market and model implied prices of European call options (on futures). Model implied prices are computed from (14), for implementing the COS method infinite summations are truncated at the 2 9 -th element and ODEs are solved numerically using an explicit Runge-Kutta (4,5) formula. 5 More precisely, we solve the following optimization problem: where is the vector of parameters,ñ K is the number of strikes,ñ T is the number of maturities, C Mkt t,k is the market price of the European call option on future while C t,k is the model price. We impose standard constraints: {V 0 , λ 0 , α, σ J , k v , θ v , σ v , k δ , σ δ , k λ , θ λ , β} > 0, −1 < ρ < 1 and k λ > β. In order to mitigate the dependence on the initial point supplied to the optimizer we randomly generate 5000 model parameters, then we evaluate the objective function in each of them and take the 10 parameters combinations with smallest objective function (the total running time for this step in our PC is around 11 minutes). Next, we run 10 different optimizations starting from those points (optimization is performed using the built in Matlab ® function fmincon with the interior point algorithm and the time required for this step is around 1 hour, i.e. nearly 6 minutes for each starting point). Finally, we take those parameters where the objective function presents the smallest value. The results of this procedure are reported in Table 1, where we show calibrated parameters, and in Fig. 4, where we compare the model and market option prices. Now, some comments are in order. The parameter controlling the mean reversion α is smaller for crude oil than the other commodities. This is consistent with literature on oil modeling. Indeed, despite mean reversion in commodity markets is a widely acknowledged stylized feature, for the specific case of oil, many authors started excluding mean reversion from price dynamics (Trolle and Schwartz, 2009;Larsson and Nossman, 2011;Shiraya and Takahashi, 2011;Cortazar et al., 2017). In particular Larsson and Nossman (2011) find no evidence of autocorrelation of log-returns from 25- May-1989to 25-May-2009. Meade (2010 finds better out of sample performances for models without mean reversion. On the other hand, we find strong mean reversion in the copper price, consistently with Schöne and Spinler (2017). For what concerns the parameters governing the volatility process, we find that the volatility of crude oil returns is much less persistent than other commodities as witnessed by the higher value of k v . Moreover, we find that the Feller condition (2k v θ v > σ 2 v ) is respected for silver and copper but not for crude oil and gold. Note that this does not indicate that the model is misspecified for crude oil and gold. Indeed, when option pricing models with CIR-type variance (e.g. Table 1 Calibrated parameters for the model in (7) Heston model) are calibrated on real option quotes, it often happens that the Feller condition is not satisfied (see e.g. Rouah, 2013, Table 6.2). The violation of the Feller condition implies the introduction of a non-negligible bias when pricing options via simulation with the Euler scheme or other similar methods (see the discussion in Begin et al., 2015). This problem is not relevant to our simulation approach outlined in Sect. 3 (see Sect. 5.3). Indeed, in our approach we simulate the variance process directly according to a non central chi-squared distribution which is always positive also when the Feller condition is not respected (see also Broadie and Kaya, 2006). The parameter ρ, which controls the leverage effect, is negative across all the commodities. This means that when the prices drop, volatility rises. The leverage effect is consistent with the phenomenon in which many investors hedge their physical risks with forward contracts. As a result, panic can break out when prices drop, pushing volatility up. Anyway, we find that this effect is much more pronounced for crude and gold than silver and copper. Regarding the convenience yield dynamics, all commodities display a similar degree of persistence. The long run mean is positive for oil and copper and negative for gold and silver. The diffusion coefficients have comparable magnitude, except for that of silver which is smaller. Finally, regarding the parameters of the jump intensity, we find that the parameter controlling the self-exciting effect β is similar among different commodities and the expected number of jumps per year (computed according to Dassios and Zhao, 2013, Proposition 2.3) is 5.98 for Crude oil, 4.35 for Gold, 6.10 for Silver and 7.26 for Copper. Finally, since there is no closed solution for A, B and D in Proposition (4), we display the values of F 1 , F 2 and F 3 together with the functions A,B,C,D, and G for a different values of u 1 ranging between 0 and 20, u 2 = u 3 = u 4 = 0 and different times τ .

Accuracy of Algorithm 1 and Asian option pricing
In order to assess accuracy of the proposed model simulation scheme we follow the procedure outlined in Cai et al. (2017, Section 4). First, we need to estimate the error we are committing Blue, red and yellow lines correspond to, respectively, τ = 1/4, τ = 1/2 and τ = 1 when simulating the model (7)-(10) using Algorithm 1. To this aim, we consider the problem of pricing an at the money (ATM) European call option on spot (for simplicity and facilitation of error comparisons, let us assume an initial price for the underlying equal for all the commodities S 0 = 100), then we compare the true option price computed from (13) using the COS method and the price of the same option obtained through a high number of simulations (following Broadie and Kaya, 2006 we use 10 9 simulations to get an accuracy at least to the fourth decimal place). The bias is the difference between the two values. Results are reported in Table 2. From this table, we can appreciate how the accuracy drastically increases with the number of time steps n. Moreover, for n = 16 we find that the bias is close to 0, meaning that if one increases further n the error will disappear. This finding also confirms that the error coming from the simulation of T 0 V s ds|V T , V 0 is negligible in practice. The speed of mean reversion α impacts on the accuracy: for crude oil we have α = 0.0637 and in this case we register the best performances of Algorithm 1, with a bias smaller (in absolute value) than 0.01 for n = 2. When α increases we need higher n to reach similar accuracy, for example, in the case of silver, we need n = 8 for an absolute bias smaller than 0.01.
Having confirmed that accuracy increases with the number of time steps n, the next step is to study the performances of Algorithm 1 in terms of trade-off between accuracy (measured through the Root mean Squared Error, RMSE = bias 2 + standard error 2 , see Li and Wu, 2019 for more details) and computational efficiency (running time expressed in seconds). Following standard literature (e.g. Broadie and Kaya, 2006, Cai et al., 2017, Li and Wu, 2019 we report the performances of the Euler scheme for benchmark comparison. First of all, it is worth noting that the proposed simulation method differs greatly from discretization methods like the Euler scheme. Indeed, we just need to split the time horizon into a very small number of time discretization steps (indeed, biases are very close to 0 for n = 16), while Euler scheme uses simple but rough approximations which work well only on small time steps, with the natural consequence that Euler scheme needs a higher number of time steps to achieve a level of accuracy similar to that of our proposed approach. However, when using Algorithm 1 and the Euler scheme for the purpose of pricing options with a limited computational budget, a trade-off is intrinsically established between increasing the number of time steps n (to reduce the bias) and the number of simulations N (to reduce the statistical error). Increasing n will improve accuracy but also increase the computational cost. We compute the bias of the Euler scheme using 10 9 simulations for different number of time discretization steps n = {200, 400, 800, 1600, 6400}. The behavior of the bias of both Algorithm 1 and the Euler scheme for different time steps is reported in log-log scale in the top panel of Fig. 6. From this figure we can appreciate how the bias of Algorithm 1 decreases faster than that of the Euler scheme. More precisely, by regressing log 10 (bias) vs log 10 (n) we get a slope around −2 (for all parameter sets) for Algorithm 1 and around −0.9 for the Euler scheme, implying that the bias of our proposed simulation method is approximately proportional to 1/(n 2 ), while the bias of the Euler scheme is approximately proportional to 1/n. Duffie and Glynn (1995) show that if the bias decays at the order of 1/(n p ), then one should increase n proportionally to N 1/(2 p) to achieve asymptotic optimality. Hence, it is possible to improve efficiency of the proposed method with a smaller computational effort on the bias reduction with respect to the Euler scheme: having estimated p ≈ 2 for our Algorithm 1 we can select n increasing at the order N 1/4 , while in the case of the Euler scheme, since p ≈ 1, n should be proportional to N 1/2 . The trade-off between accuracy and computational efficiency in a log-log scale is shown in the bottom panel of Fig. 6 while numerical values of RMSE and running times are reported (along with biases) in Table 3. Results show that the RMSE of the proposed simulation method decays faster than the Euler scheme. In particular, we get the following convergence rate of the RMSE for Algorithm 1: 0.44 for crude oil, 0.43 for gold, 0.44 for silver and 0.42 for copper. These are only slightly smaller than those of a theoretically unbiased estimator (which would be around 0.5, see e.g. Broadie and Kaya, 2006) but higher than those of the Euler scheme, which are around 0.30 for all the parameter sets considered, confirming the superior performances of the proposed approach with respect to the benchmark. In Fig. 6, in the case of silver we also included the case with n = 8 and N = 1024 × 10 4 for Algorithm 1 to better highlight its superior performances. Another   Table 1. Other parameters: S 0 = 100, K = 100 and T = 1. Further notes: see Table 2 important aspect in the comparison between the proposed approach and the Euler scheme is that in the latter the discretization of the variance process may generate negative values in intermediate steps with a significant probability when the Feller condition 2k v θ v > σ 2 v is violated (this is the case of the calibrated parameters for crude oil and gold). This brings extra error in the simulation procedure, explaining the poor performances of the Euler scheme in the case of gold (where, in addition, the initial value of the variance process V 0 is very small, increasing the probability of the variance process reaching zero). Finally, given model parameters and having tested the accuracy of the simulation scheme, we consider the problem of pricing an Asian option using Algorithm 2. We assume an initial price S 0 = 100 and compute the price of European and Asian call options (geometric and arithmetic averaging) using (13) and (23) with discrete monitoring (12 monitoring dates) and several strikes. ODEs are solved numerically using an explicit Runge-Kutta (4,5) formula, while infinite summations for the implementation of the Fourier-Cosine method are truncated at the 2 10 element. Algorithm 2 is implemented with N = 10 6 and N CV = 10 4 simulations. Results are reported in Table 4. European call option price is computed in approximately 1 second, the geometric Asian option price is computed in around 8 seconds, while the price of the arithmetic counterpart is obtained in around 110 seconds for all the parameter sets. The usage of the geometric Asian option price as control variable turns out to be extremely useful, allowing for a variance reduction around 99% across all the parameter sets, maturities and strikes.

Conclusions
In this paper we propose a new model for pricing commodity options. The model accounts for mean reversion, stochastic convenience yield, stochastic volatility and stochastic jump intensity. For the latter, we provide empirical evidence of self-excitation across four different commodity markets. After presenting the model under the historical measure, we introduce a structure preserving change of measure and describe the model under a risk-neutral measure. Then, by calibrating the proposed model on real market option prices, we find an excellent fit. We develop an efficient simulation scheme for the proposed model, we identify sources of error and present a comparison with the classic Euler scheme. Finally, we derive semi-closed formulas for the price of geometric Asian options under the proposed model and combine this result with the simulation scheme to develop a Control Variate simulation strategy for accurate evaluation of arithmetic Asian option prices, which are very popular derivative instruments in commodity markets. The methodology is able to provide accurate results at whereJ x dN t denotes the non-compensated jump process. By this way, (A.6) becomes independent of the others. The last equation, describing the volatility dynamics, is the classical equation of the CIR type, for which existence and uniqueness of a strong solution can be proved by applying Yamada-Watanabe theorem (see Karatzas and Shreve, 1991, Proposition 2.13, Chapter 5, p. 291). As a continuous process, on every interval [0, T ] it will be bounded. (A.3) is an Ornstein-Uhlenbeck stochastic differential equation driven by a Wiener process and existence and uniqueness of the solution of this equation is a classical result of stochastic calculus. As a continuous process, on every interval [0, T ] it will be bounded as well.
Finally, we can focus on (A.5), describing the log-returns dynamics. By performing an equivalent change of measure with the following likelihood process: the jump term appearing in (A.5) turns into a Compound Poisson process with constant intensity. We point out that this measure change involves only the jump term, since all the other terms are assumed to be independent on the jump part. Moreover the measure change will establish a one-to-one mapping between the processes under the two measures. The likelihood process L t , under the assumption E P [e J x ] can be proved to be a true martingale since E P [ t 0 λ s ds] < +∞, ∀t ∈ [0, T ]. (A.5) is now a linear stochastic differential equation driven by a Lévy process, and the classical result of existence and uniqueness of strong solutions holds (a proof can be found in Applebaum (2004, Theorem 6.2.3, p. 304), since the boundedness of both δ t and V t on the time interval [0, T ] ensures the required Lipschitz conditions. with terminal conditions A(T , T ) = 0, B(T , T ) = w, C(T , T ) = 1, and denote by = e κ 1 (ξ )β ψ(ξ ), we can write: By inserting these expressions into (C.2) we obtain the following equation, that must be satisfied for every value of λ, m: We obtain C(t, T ) = 1 and Finally, by using conditions 6, we get: We can repeat the same computation of the conditional moment generating function of λ under P without the introduction of the terms m t and m T and we get the following system of ordinary differential equations: By comparing the equations written under Q and P, the result follows. We can perform the computations of the conditional moment generating functions for all the remaining state variables X t , V t , δ t (or alternatively of the joint conditional moment generating function), by following the procedure outlined above for λ t , both under Q and P, and, by direct comparison of the corresponding equations, the relations among the parameters under Q and P can be obtained. We omit all these computations, which are tedious and straightforward. In Appendix D an explicit solution is computed for the joint conditional moment generating function under Q.
By separating the state variables we arrive at the ODE system of the thesis.